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ABSTRACT 

We study the interaction of a proto-planetary disk and a planet on a highly inclined orbit in the 
linear regime. The evolution of the planet is dominated by dynamical friction for planet masses 
above several Earth-masses. Smaller planets are dominated by aerodynamic drag, especially 
for very high inclinations and retrograde orbits. 

The time-scales associated with migration and inclination damping are calculated. For 
certain values of the inclination, the inclination damping time-scale is longer than the migra- 
tion time-scale and the disk lifetime. This result shows that highly inclined planets can not 
(re-)align with the proto-planetary disk. 

We discuss the dependence of numerical simulations on the gravitational softening pa- 
rameter. We find only a logarithmic dependence, making global three dimensional simulations 
of this process computationally feasible. 

A large fraction of Hot Jupiters is on highly inclined orbits with respect to the rotation axis 
of the star. On the other hand small-mass planetary systems discovered by the Kepler mission 
have low mutual inclinations. This shows that there are two distinct formation mechanisms at 
work. The process that creates inclined Hot Jupiters does not operate on small mass planets 
because the damping timescales are so long that these systems would still be inclined today. 

Key words: Solar System; planets and satellites: formation; celestial mechanics, accretion, 
accretion discs, methods: analytical, methods: numerical 



1 INTRODUCTION 

The number of confirmed extra-solar planets was at 556 when this 
paper was submitted 1 . Many of them, mostly the high-mass, close- 
in Hot Jupiters, are on rather extreme orbits, especially when com- 
pared to the solar system. Using the Rossiter-McLaughlin effect, 
several planets have been found to be in highly inclined orbits with 
respect to the sky-projected spin axis of the star (Triaud et al. 2010; 
Simpson et al. 201 1, and references therin). Although still debated 
(Lai et al. 2011), it is generally assumed that the proto-planetary 
disk is roughly aligned with the spin axis of the host star. 

Several mechanisms have been proposed to account for plan- 
ets on orbits with high inclinations (e.g. Fabrycky & Tremaine 
2007; Chatterjee et al. 2008). Here, we do not address the question 
on how the planet got on the orbit in the first place. 

Before highly inclined planets had been found, one of the 
main areas of theoretical work was the interaction of planets and 
proto-planetary disks. The robust result is planetary migration, al- 
though the precise speed and also the direction are still being de- 
bated and most likely depend strongly on the global disk structure 
(Paardekooper et al. 2010, 2011). This naturally leads to the ques- 
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1 By the time this paper was accepted the number has gone up to 759. 



tion of how a planet on a highly inclined orbit interacts with a proto- 
planetary disk. Several authors have looked at the evolution of mod- 
erately inclined planets in three dimensional simulations (Marzari 
& Nelson 2009; Bitsch & Kley 201 1 ; Cresswell et al. 2007) 

In this paper, we consider a planet that is on a highly inclined 
orbit, not embedded in the disk, and crossing a vertically stratified 
disk twice per orbit. The planet creates a density perturbation in the 
disk. These perturbations interact gravitationally with the planet. 
This interaction is known as dynamical friction. Furthermore, aero- 
dynamical drag from accreting material might become important. 

In section 2 we first calculate the linear response of the disk 
and the resulting friction force on the planet analytically. We then 
estimate the relevant time-scales and compare the effect to aerody- 
namical drag. In section 3 we describe numerical simulations ver- 
ifying the analytic result. We finally discuss the consequences in 
section 4. 



2 LINEAR THEORY 

We assume that the disk is in Keplerian rotation with angular ve- 
locity Q. at a distance a from the host star with mass M* . We further 
assume that the planet is on a circular orbit with velocity v k = On 
and inclination ;'. The relative velocity between the disk and the 
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i/2 



Figure 1. Illustration of the setup; a planet with inclination i passing 
through the disk on a highly inclined orbit from top left to right bottom. 



planet is then v imp = 2v k sin(7/2). The force due to dynamical fric- 
tion or aerodynamic drag felt by the planet is anti-parallel to the 
velocity of the planet relative to the velocity of the disk. This setup 
is illustrated in Fig. 1. The force has an angle i/2 to both the normal 
direction of the planet's orbit and the normal direction of the disk. 
We ignore the shear and curvature of the disk, and assume that there 
is no net force in the radial direction. Without loss of generality, we 
will work in units where v k = 1 , a = 1 , GM* = 1 unless otherwise 
stated. 



2.1 Dynamical friction 

When the planet comes close to the disk plane, it induces density 
perturbations in the disk. These then exert a gravitational force on 
the planet. The force is directed so that the relative velocity of the 
disk and the planet always decreases. Thus the name, dynamical 
friction. 



2.1.1 Linearized equation 

In analogy to Ostriker (1999) we start by defining the perturbed 
density a and velocity fi by p(x, t) = po(x)[l + a(x, t)] and v(x, t) = 
c,/?(x, t). We consider an isothermal equation of state with constant 
sound-speed c s . The linearized continuity and momentum equa- 
tions in a frame moving with the disk are given by 



da 



+ c s V/3 = 



dp 
dt 



+ c,Vff + c/VO = 



(1) 

(2) 



where O = -Gm/\x - x p \ is the potential due to the planet 
at position x p . The linear equations are valid in the limit 
that ar,|/3|,|Vpo/pol « 1. As shown by Kim & Kim (2009), this 
is justified as long as the dimensionless nonlinearity parameter 77 is 
less than unity. For an Earth-like planet on a 90° orbit we have 

Gm 

n = 



c 2 r p (M 2 - 1) 



~ 0.03. 



(3) 



The background steady state is a stratified disk in the z direction 
with po(x) = p exp[-z 2 / H 2 ], where H is the thickness of the disk. 



In a Keplerian accretion disk we have H = V2 c s Cl 1 . By eliminat- 
ing p, the above equations lead to 

1 d 2 a 



V 2 a- 



dt 2 



--^v 2 o. 

c% 



(4) 



This can be solved with a Green's function (Ostriker 1999). The 
steady state solution is given by 



a(x, t) = 



Gm/c 2 



^s 2 +R\\-M 2 ) o otherwise . 



1 M< 1 

2 M > I and s/R < - VM 2 - 1 



The above equation describes a Mach cone in the supersonic case 
with no perturbations outside of the cone (see Fig. 1). s is the co- 
ordinate in the direction of the Mach cone and R is in the direction 
perpundicular to it. Because these equations are linear, we can eas- 
ily generalize this result to a softened planet potential of the form 
(£>(R, s) = -Gm/ y/R 2 + s 2 + b 2 with softening length b. For this 
potential, the corresponding mass distribution is given by 



P P {R, s) = 



3mb 2 



An (b 2 +R 2 + s 2 f 2 ' 
Thus, the density perturbation created by this mass distribution is 

Pp(|x'|) , , 



(5) 



a(x, t) 



Iff 



a(x - x , t) 



(6) 



As we will integrate over x in the following section, we can simply 
use Fubini's theorem and do not need to evaluate the integral here. 



2.1.2 Change in orbital energy after one passage through the 
disk 

The specific gravitational force exerted on the planet from an un- 
perturbed, infinitely large stratified sheet with density profile p„(x) 
is symmetric to the mid-plane. Because of this symmetry, there is 
no net change in orbital energy of the planet after one crossing. 

However, the force from the perturbed density a is not sym- 
metric, the planet therefore changes its specific orbital energy af- 
ter each passage trough the disk. This can be calculated as an 
integral over the specific force F along the path of the planet 
x„ = (x P ,y p ,z p T l , 



AE 



I 



= F(x„(0)< 





cos(j') 
sin(i) 



dt. 



Using the definitions from the last section (see also Fig. 1), we can 
change the coordinate system and describe any point x = (x,y, z)~' 
with R, s and </> using the relations 

x = x p + R sin(0) 

y = y p + R cos(0) cos(('/2) + J- sin(i'/2) 
z = z p + R cos(0) sin(;'/2) + s cos(i/2). 
The force on the planet is then given by 
p(x) (x„ - X) 



^-ojff 

-Iff 



d'x 



\(x p -x) 2 + b 2 \ 3 ' 2 
2G 2 mp /c 2 s __ 2/H 2 Xp-x 
V-v 2 +tf 2 (l-M 2 ) \ s 2 +R 2 + b 2 



3/2 



R dR ds dep. 
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The change in orbital energy can then be integrated over <f> and the 
trajectory of the planet, such that 



-2 --- /H- 

mp c/e -'»' 



V-S' 2 +tf 2 (l-M 2 ) |. V 2 + R 2 + b 2 



3/2 



R dR ds d(p dx p 



= - 4n 3 ' 2 HG 2 mp /c 2 A(M, r>, b)/M 2 , 

' sin(;j 

where we have defined the remaining two integrals as 

M 2 s R 

A(M,r, ' 



J J \ S 2 + 



[s 2 +R 2 (l-M 2 )] 1/2 \s 2 + R 2 +b 2 



3/2 



dRds. 



The integral has the form of a Coulomb logarithm and we have 
introduced an additional cut-off scale at large distances (e.g. the 
disk size), r > . This allows us to solve the Coulomb integral exactly, 

A(M, r>, b) = - ^M 2 - 1 • atan T> 

i V(M 2 - 1)(* 2 + rl) 



log 



(7) 



Note that if we further introduce the surface density S = H V^Po 
and substitute the expression for the Mach number, we can rewrite 
the specific change in orbital energy as 



AE„ = — 



nG 2 mL 



v: sin(;'/2) sin(i') 



A(M, r> ,b). 



(8) 



The result is almost independent of H , c s and M. The only depen- 
dencies are hidden in the function A which behaves like a Coulomb 
logarithm. A is typically of order — 1 — 10 and reflects the uncer- 
tainty about the small scale processes near the planet's surface and 
the large scale effects at scales where the disk curvature and shear 
become important. Equation 8 reduces to that of dynamical friction 
in a collision-less medium in the limit of r> » b and M » 1 . To 
see this note that we can take that limit in Eq. (7) to get 



A(M, r > ,b) — > log 



(9) 



which should be compared to Eq. (8.5) in Binney & Tremaine 
(2008). 

It is also worth noting that our result is qualitatively similar 
to those summarized by Artymowicz (1994) (see his Eq. 69) but 
includes explicitly the dependence on the gravitational softening 
parameter b. Thus, it can be used as a benchmark in numerical sim- 
ulations that make use of such a softening. Furthermore, as we will 
show later, we can argue that a softening parameter that is only 
slightly smaller than the disk scale height can be used in numerical 
simulations with a moderate error. 



2.1.3 Time-scales 

The rates of change in semi-major axis and inclination can be ex- 
pressed in terms of the specific normal A' and tangential force f 
(Burns 1976). For a circular orbit, these are 



da 
dt 



= 2 



GM+ 



a T 



and 



di 
dt 



GM t 



N. 



Averaging over one orbit, using Eq. (8) and noting that the force 
from dynamical friction is inclined by an angle i/2 with respect to 
the normal of the orbit, we get 

4nmLa i 1 



M\ sin(i72) sin(0 



- A(M, r> , b), 



4a 2 

Aa= — AE a 

2a « , 4mnLa 2 1 

Ai = AE a tan- 1 0'/2) = = A(M, r>, b). 

GM„ Ml 4sin 3 (;72) 

The time-scales for inclination and semi-major axis damping are 
then 

M 2 

r„ = sin(i/2)sin(() A-\M, r> ,b) Or 1 , (10) 



2ml.a 2 
M 2 

T t = Vsin 3 (i'/2) A _1 (M, r > ,b) Or 1 . 

mLa 2 



(11) 



2.2 Aerodynamic drag 

The aerodynamic drag due to mass accretion might also be impor- 
tant. The accreted mass after one crossing is 



Am = 



t» 2 



(12) 



cos(i/2) 

The accretion is dominated by the geometrical cross section as long 
as the Bondi radius r a - 2Gm/v 2 mp is small compared to the physi- 
cal radius of the planet r p . For a 100 Earth mass planet with a size 
of 10 Earth radii on an orbit with ( = 90° we have r a ~ 0.7r p . Thus 
the geometrical term is dominant. We will therefore assume this to 
be the only term for the remainder of this paper. The direction of 
the resulting force is the same as for the dynamical friction, i.e. 45° 
in the z<p plane for a planet on a 90° orbit, and the magnitude of the 
specific force is 

Q. Am 

fd = V imp . 



(13) 



The relevant time-scales are then 



- c l - _L l GM * 
" a 2T y a 

i i IGM+ 



(14) 



(15) 



where the tangential and normal components of the force can be 
decomposed as in the previous section. Putting everything together, 
we have 

m cos(//2) 



nirj 2tan(i'/2) 

111 

?cos(//2). 



n» 2 



(16) 



(17) 



Note that the dependencies on mass and inclination are different 
compared to the dynamical friction time-scales. 

2.3 Comparison of time-scales 

As an illustration, let us assume a solar nebula with a surface den- 
sity S = 4200 g/cm 2 (Weidenschilling 1977) and a 10 Earth-mass 
planet at 1 AU. We will also set A = 4 for simplicity. The time- 
scales given by Eqs. (10) and (1 1) as well as the time-scales for the 
aerodynamic drag, Eqs. (16) and (17), are plotted in the first panel 
of Fig. 2 as a function of the inclination i. We also plot the time- 
scales for type I migration calculated by Tanaka et al. (2002) and 
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2 100000 - 



semi-major axis (dynamical friction) - 

nclination (dynamical friction) 
semi-major axis (aerodynamic drag) 
inclination (aerodynamic drag) 
semi-major axis (Tanaka) 
nclination (Tanaka) - 



inclination [degree] 

(a) 10 Earth mass planet. 




-major axis (dynamical friction) 

nclination (dynamical friction) 

major axis (aerodynamic drag) 

inclination (aerodynamic drag) 
semi-major axis (Tanaka) 
nclination (Tanaka) 



100 



120 140 



(b) 100 Earth mass planet. 
Figure 2. Migration and inclination damping time-scales. 



Tanaka & Ward (2004) as a comparison. Their inclination damping 
time-scale is only valid for small mass planets embedded in the disk 
(i < few degrees). Bitsch & Kley (201 1) showed that this is roughly 
consistent with numerical simulations when planets are completely 
embedded in the disk. 

In the lower panel of Fig. 2, we plot the same quantities for a 
100 Earth mass planet. Here, the aerodynamic drag is reduced by a 
factor VlO, whereas dynamical friction is increased by a factor 10 
compared to the 10 Earth mass case. 

For a 10 Earth mass planet which is not embedded in the disk, 
the dynamical friction dominates the evolution of the semi-major 
axis up to about 50° inclination. The evolution of the inclination is 
dominated by dynamical friction up to about i ~ 70°. Beyond that, 
the aerodynamic drag becomes important. 

For a 100 Earth mass planet, dynamical friction dominates up 
to retrograde orbits with i ~ 130°. For orbits above i ~ 100°, the 
inclination damping time-scale is longer than the migration time- 
scale. This leads to interesting questions about the orbital evolution. 
From this analysis alone, we would expect an excess of planets with 
several tens of Earth masses on highly inclined orbits. 



3 NUMERICAL SIMULATIONS 

We ran high resolution, three dimensional hydro-dynamical simu- 
lations of a planet interacting with a proto-stellar disk to confirm 
the validity of the linear theory described above. In this paper we 
only simulate a small wedge of the disk to speed up calculations 




(a) Tangential force on the planet after removing the contribution due to the 
unperturbed disk as a function of time for three different inclinations. The 
mass of the planet is one Earth mass. 




(b) Fractional change in the planet's orbital energy after one crossing as a 
function of the inclination. The mass of the planet is one Earth mass. 




softening lengtO b [AU] 

(c) Fractional change in the planet's orbital energy after one crossing from 
numerical simulations and linear theory as a function of the softening length 

Figure 4. Results from the hydrodynamic simulations. 



as we are only interested in one passage of the planet through the 
disk. The results presented in this section should be seen as a proof 
of concept, illustrating that it is possible to have a well defined 
and converged numerical simulation. We will study interesting non- 
linear effects such as the question if the planet is able to open a 
gap and where the precise transition between the Tanaka & Ward 
(2004) regime and dynamical friction is, in much more detail in an 
upcoming paper. 

We use the Athena code in cylindrical coordinates with a res- 
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0.975 0.9S 0.985 0.99 0.995 1 1.005 1.01 1.015 1.02 1.025 



0.975 0.98 0.985 0.99 0.995 1 1.005 1.01 1.015 1.02 1.025 



Figure 3. Comparison of linear theory (left) and numerical simulations (right). Plotted are slices of the normalized density p(x)/po(x) = a(x) + 1 for an Earth 
mass planet crossing the disk plane. The snapshot is taken when the planet crosses the mid-plane of the disk (z = 0) and along the direction of the force (see 
Fig. 1). 



olution of up to (N r ,Nj,N-) = (768, 1024,256). The size of the 
wedge is typically (L r , L$, L z ) = (0.6, 0.8, 0.2). An isothermal equa- 
tion of state is used with a sound-speed of c s = 0.05. The planet is 
on a i = 90° orbit so that the Mach number is M = 28.3. Note 
that because of the CFL condition and M » 1 the time-step in the 
simulation is very small, typically ~ 10~ 5 . 

Fig. 3 shows a slice of the density perturbation p(x)/p (x) = 
o-(x) + 1 in the direction of the force (45° in the z<p plane for a 
planet on i = 90°). This is given analytically by Eq. (6) on the left 
and by the numerical simulation on the right. The overall agreement 
between the linear calculation and the numerical result is very good 
with two differences worth noting. First, the numerical solution is 
not symmetric with respect to r = 1. This is due to the curvature 
of the orbit, whereas the analytic model on the left assumes that 
the orbit is on a straight line 2 . Second, the numerical solution is 
slightly smeared out compared to the analytic solution. This is due 
to numerical viscosity. 

We show the outcome of simulations with a one Earth mass 
planet on an orbit with ( = 45°, 90°, 155° in Figs. 4 (a) and (b). 
Fig. 4 (a) shows the tangential force onto the planet after removing 
the part from the unperturbed disk as a function of time. The planet 
is crossing the mid-plane at / = 0.25. The differences in the force 
can be understood as follows. The interaction of the planet with the 
disk in the i = 45° case is long and strong (low M). For the i = 90° 
case the interaction is shorter and weaker. For the i = 155° case the 
interaction gets longer again (because the planet spends more time 
in the disk) but even weaker (high M). 

Integrating this component of the force along the trajectory 
of the planet gives the change in orbital energy, as calculated in 



Note that the plot is stretched in the z<p direction. 



Eq. (8). We measure the change in specific orbital energy of the 
planet in each of those simulations. The results are plotted in 
Fig. 4 (b), together with the theoretical predictions. In the / = 45° 
case there is a noticable difference which is most likely due to the 
curvature of the disk that becomes more important for small i. 

Finally, let us focus on the effect of the softening length. The 
relevant scale in the problem is H 1 la ~ 0.0025. However, the ex- 
pression for A in Eq. (9) suggests that there is only a weak de- 
pendence on the softening parameter b. We therefore ran multiple 
simulations with different softening parameters to verify this. In 
Fig. 4 (c) we plot the fractional change in orbital energy of the 
planet after one disk crossing as a function of the softening param- 
eter. The solid line is given by Eq. (8). The agreement is very good, 
over more than one order of magnitude. Effects that haven't been 
considered in the linear calculation, such as shear and curvature of 
both the orbit and the disk, are all 0(H) ~ 0.05 or smaller, leading 
to the small error seen in Fig. 4 (c). 



4 DISCUSSION 

In this paper, we derived the linear equations for the evolution of 
highly inclined planets interacting with a proto-planetary disk. De- 
pending on the mass of the planet and the inclination, either dy- 
namical friction or aerodynamic drag are dominant. For certain in- 
clinations, the migration time-scale is shorter than the inclination 
damping time-scale and even the disk lifetime. 

For very small inclinations (i < 5°), when the planet is embed- 
ded in the disk, other effects become important. The time-scales 
given by dynamical friction approach values comparable to stan- 
dard type I migration as estimated by Tanaka & Ward (2004). In 
reality there is a smooth transition between the two regimes (see 
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Fig. 2), although this is not captured by the calculation presented 
above. 

Using a completely different approach, Papaloizou & Lar- 
wood (2000) calculate the damping timescales of highly eccentric 
planets. Their eccentricity timescale scales as ~ e 3 . We have a scal- 
ing of ~ i 4 for small i (see Eq. 11). We can understand the difference 
by noting the interaction for an inclined planet can only take place 
while it is embedded in the disk, whereas a highly eccentric planet 
is always embedded in the disk. 

A large fraction of high-mass, close-in planets (Hot Jupiters) 
has highly inclined orbits (Triaud et al. 2010). On the other hand, 
low mass planet candidates in multiplanetary systems discovered 
by the Kepler mission seem to be aligned with low to moderate mu- 
tual inclinations (Lissauer et al. 2011). If this trend is correct, then 
we can use the results from our calculation to predict that there are 
two distinct formation scenarios for low and high mass (and close- 
in) planets. The timescale for realignment of small mass planets is 
very long. Therefore these objects have never been inclined other- 
wise we would still see them on inclined orbits today. Whatever 
process moved Hot Jupiters on inclined orbits does not operate on 
small mass planetary systems. 

This calculation might also have interesting implications for 
other systems where dynamical friction occurs in the presence of 
a gaseous disk, such as for binary stars and super-massive black 
holes (e.g. Kim 2010). 

In a follow-up paper we will study the long-term and non- 
linear evolution with detailed global hydro-dynamical simulations. 
Using the results of this paper, a softening length of just slightly 
smaller than the disk scale height is sufficient to capture the relevant 
physical process. Without these results, numerical simulations of 
planets on highly inclined orbits interacting with a proto-stellar disk 
would not be feasible. 
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